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ABSTRACT 

This work provides an in-depth mathematical description of the response functions 
that are used for spatial and spectral analysis of X-ray data. The use of such functions 
is well-known to anyone familiar with the analysis of X-ray data where they may be 
identified with the quantities contained in the Ancillary Response File (ARF), the Re- 
distribution Matrix File (RMF), and the Exposure Map. Starting from first-principles, 
explicit mathematical expressions for these functions, for both imaging and dispersive 
modes, are arrived at in terms of the underlying instrumental characteristics of the tele- 
scope including the effects of pointing motion. The response functions are presented in 
the context of integral equations relating the expected detector count rate to the source 
spectrum incident upon the telescope. Their application to the analysis of several source 
distributions is considered. These include multiple, possibly overlapping, and spectrally 
distinct point sources, as well as extended sources. Assumptions and limitations behind 
the usage of these functions, as well as their practical computation are addressed. 

Subject headings: methods: data analysis, analytical — X-rays: spectra 

1. Introduction 



It is a basic tenet of X-ray spectral analysis that the source flux incident at the telescope 
is related to the observed count rate through an integral equation involving the effective area of 
the telescope. The most commonly accepted technique for dealing with this equation involves 
the use of spectral analysis programs such as xspec (Arnaud 1996). The effective area is input 
into these programs via a file called the Ancillary Response File, or ARF. In addition, the energy 
resolution of the detector is specified by the Redistribution Matrix File, or RMF. This work presents 
formal descriptions of the quantities embodied by the ARF and RMF in terms of the underlying 
instrumental responses, making a clear connection between the incident source flux and the observed 
count rate. 

Roughly speaking, the effective area of an X-ray telescope composed of a mirror and a detector 
is more or less the product of the effective area of the mirror with the quantum efficiency (QE) 



- 2 - 



of the detector. However, the ARF, which relates observed counts to a source flux, also depends 
upon observation-dependent quantities such as the detailed aspect history of the telescope, its point 
spread function (PSF), and upon details of the analysis itself, e.g., the filtering and binning of the 
observed data. 

With the advent of the Chandra X-ray Observatory (CXC 2000), all of the subtleties intro- 
duced by the shape of the PSF and the telescope aspect motion were deemed important for the 
computation of an ARF. Chandra is well calibrated (Weisskopf et al. 2000) and, with its unprece- 
dented combined spectral and spatial response, a precise definition of the ARF that incorporates 
the effects of spacecraft motion and the PSF is necessary in order to perform spectral analysis at 
the highest resolution of the instrument. 

The main goal of this work is to present an explicit first-principles derivation of the ARF 
by including the proper treatment of telescope motion (e.g., dither) and the PSF, as well as any 
filtering of the data. The bottom-up approach taken here necessarily implies that a meaningful and 
consistent derivation can be achieved by considering the role of the ARF in spectral analysis. As 
a result, the ARF is presented in the context of integral equations that connect the incident X-ray 
flux to an expected count rate. 

Traditionally, the ARF and RMF have been used primarily for the analysis of spectral image 
data. An important aspect of this work is to extend this approach to the analysis of dispersive 
spectral data, such as data obtained by Chandra or Newton. To this end, definitions of an ARF 
and an RMF are presented that are suitable for dispersive data analysis, and which may be utilized 
by existing spectral analysis software. 

One of the original motivations for this work was the need to create a related object, an 
exposure map, for use in the analysis of data obtained by the Chandra X-ray Observatory. This 
paper also gives a rigorous definition of the exposure map and discusses some of its uses and 
its limitations in spectral image analysis. The resulting definition is consistent with current use, 
intuition, and physics. 

The next section contains a discussion of the general response of an X-ray telescope and also 
serves to introduce the notation and conventions used throughout this paper. Although originally 
inspired by the need to create ARFs and exposure maps for Chandra, the presentation has been 
kept as general as possible without focusing on any particular telescope or instrument. A derivation 
of the imaging ARF follows in section 3 where its application to several problems is considered. 
These include the problem of multiple overlapping point sources. Section 4 contains a definition of 
the exposure map and explores its use as well as its limitations in dealing with extended sources. 
The definition of a dispersive ARF and RMF that are suitable for use in the analysis of dispersive 
spectral data are given in section 5. Following the summary of the paper is an appendix that 
considers the practical computation of these objects. 
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2. Telescope Response 

Let S(X,p)d\ dp represent the number of photons per unit area per unit time incident upon 
the telescope with directions that lie in the cone of directions between p and p + dp, and whose 
wavelengths lie between A and A + dX. This source flux is assumed to be time-independent 1 . 
Similarly, let So(h,a,t)da represent the expected number of counts per unit time, in pulse-height 
bin h, within a region of area da at the position a on the detector. (Although in this paper, h 
is a discrete quantity that represents a pulse-height channel, one could easily generalize ft to a 
continuous quantity such as a voltage by introducing the infinitesimal dh and replacing J2h by 
fdh.) 

In general, the count rate will be time-dependent even if the source flux does not vary with 
time. This time-dependence arises from several effects, including but not limited to telescope 
pointing motion (e.g., pointing wobble or dither), detector electronics, telemetry saturation, and 
thermal expansion effects that cause the movement of individual telescope subsystems with respect 
to one another. 

As used here, pointing is the measured, generally periodic movement of a coordinate system 
attached to the optical axis with respect to a coordinate system fixed with respect to the sky. These 
two coordinate systems are related by a time-dependent rotation matrix TZ(t) which completely 
characterizes the pointing. It is assumed that complete knowledge of lZ(t) is available from, e.g., 
an on-board aspect system, at some level of accuracy (see section 2). 

In the coordinate system fixed with respect to optical axis, the source flux will appear to be 
time-dependent according to 

S(X,p,t) = S(X,K- 1 (t)-p). (1) 

This equation simply expresses the fact that an observer fixed to the telescope will see a time 
dependent source and that this induced time dependence is a direct consequence of telescope motion. 
In the telescope coordinate system, the source flux and the count rate are related to one another 
via the equation 

S D (h,a,t) = J dxj dpR T (h,X,a,p,t)S(X,K-\t)-p), (2) 

which defines the total response Rr(h, X,a,p,t) of the telescope. It is an extremely complicated 
function that incorporates all elements of the telescope such as the detector and its electronics, 
the mirror, and diffraction gratings, if present. Since the units of So(h,a,t) are counts per unit 
detector area per unit time, and the units of S(X,p)dX are photons per unit aperture area per unit 
time, it follows that the response Rr(h, X,a,p,t) is a unitless quantity (counts/photon x aperture 
area/detector area). 



1 If the source flux varies in time such that its time-dependence is uncorrelated with the spatial and spectral shape, 
as is often the case, then one can always factor out the time-dependence and handle it separately. This technique 
is discussed in more detail below. The treatment of more complex time-dependent sources that do not admit this 
factorization is beyond the scope of this paper. 
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The mathematical form of equation (2) represents a linear mapping from S to So- Hence, 
strictly speaking, equation (2) is valid only when one can neglect non-linear effects such as local 
gain depression or pile- up, i.e., a non-linear effect caused by the finite time resolution of the detector. 
Pile-up can be treated, at least in principle, by making the response Rr(h, X,a,p,t) a function of 
the incident flux S(X,p). However, this topic is beyond the scope of the present work and will be 
addressed elsewhere (Davis 2000). 

The explicit time-dependence of Rr(h, X,a,p,t) is due to effects associated with telemetry 
saturation, thermal expansion, and so on. It is assumed that all but the time-dependence associated 
with internal movement of the telescope subsystems may be encapsulated in a function T(a, t) that 
factors out of Rr(h, X,a,p,t), i.e., 

R T (h,a,X,p,t) = T(a,t)R(h,a,X,p,t). (3) 

In this equation the residual time-dependence of the response function R(h, a, A, p, t) depends only 
upon the relative movements of the individual telescope subsystems. The function T(a,t) can 
be thought of as representing the the so-called "good-time intervals", or GTIs. However, T(cr,t) 
could play a more general role than this because it could include the effects of bad-pixels 2 , which 
themselves may be time-dependent, as well as any dead-time factors associated with telemetry 
saturation. It should also be noted that although the incident source flux has been assumed to be 
independent of time, T(a, t) could encompass time-dependent variations in the source flux as long 
as the spectral shape itself does not depend upon time. In this case, only the amplitude of the flux 
is time-dependent and this time-dependence may be factored out of the incident source flux and 
into the function T(a,t). Further exploration of this possibility is left to the reader. Suffice it to 
say that the actual values that T(a, t) take on are not important in what follows. Hence, equation 
(2) can be written 

S D (h,a,t) = T(a,t) J dX J dp R(h,a, X,p,t)S(X,TZ' 1 (t) ■ p) (4) 

= T(a,t) J d) ^J dp R(h,a,X,K(t) ■ p,t)S(X,p), 

where the last form follows from the change of variables p — ► Tt{t) • P and noting that because the 
matrix TZ(t) is orthogonal, the Jacobian of the transformation is unity. 

Assuming that the response depends upon several independent telescope subsystems, it may 
be factored it into subsystem-dependent pieces. The specific form of the factorization will depend 
upon the actual physical relationships between the various subsystems. For a prototypical X- 
ray telescope with a focusing mirror at the aperture of the telescope that focuses X-rays onto a 
position-sensitive detector, an appropriate factorization is given by 

R(h,a,X,p,t) = D(h,a,X) J dp' 5(a - a(p', t))T(\, p',p)M(X,p). (5) 



2 As described below, the effects of static bad-pixels are assumed to be contained in the detector response function 
D(h,a, A). 
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A different factorization may have to be use to describe the response of some other type of telescope 
(see section 5 for the factorization needed to to describe the presence of a diffraction grating). In 
this equation, M(A, ft) is the off-axis effective area of the mirror, and D(h,a,X) represents the 
probability that a photon with wavelength A at position a on the detector will give rise to a pulse- 
height h. The term involving the delta function symbolizes the passage of a photon with direction 
ft' from the mirror to the position a on the detector via the coordinate transformation a(ft',t). In 
general, this function is time-dependent and represents any relative motion that exists between the 
mirror and the detector. 

The PSF (Point Spread Function) of the telescope is represented by the function T(\,p' ,ft), 
which is assumed to satisy the normalization condition 

1 = J dft'F(\,ft',ft). (6) 

Its definition is based upon the idea (see figure 1) that the mirror itself (e.g., one of a typical type-I 
Wolter design) may be modeled by the appropriate probability distribution for photons to enter 
and leave the mirror at a single point (the so-called "mirror node"), with an effective area given by 
M(X,ft). As indicated in figure 1, the PSF defined in this way will also depend (implicitly) upon 
the position of the detector and its relationship to the focal surface. For simplicity, in the following 
it will be assumed that the variations in the movement of the detector with respect to the mirror 
during the course of an observation are small enough that the PSF may be regarded as independent 
of time to sufficient accuracy. The reader should note that a perfect PSF defined in this sense is 
represented by 

•^perfect (A, ft', ft) = 5(ft' - ft). (7) 



Generally speaking, the detector response D(h,a,X) can be factored into a QE (Quantum 
Efficiency) Q(a, A) and a redistribution function DrQi, a, A): 

D(h,a,X)=D R (h,a,X)Q(a,X). (8) 

The function D R (h,a,X) is known as the redistribution matrix function, or RMF (George et al. 
2000), and represents a mapping, or redistribution, from wavelength to pulse- height by the detector. 
Without any loss of generality, it is assumed to be normalized to unity via 

l = Y J D R {h,a,X). (9) 

h 

Alternatively, the quantum efficiency may be defined by 

Q(a,X) = Y,D(h,a,X). (10) 

h 

In general, as indicated here the RMF varies with position on the detector, although many appli- 
cations assume a spatially constant RMF. In contrast, the quantum efficiency function Q(a,X) is 
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(a) 



Detector 




Fig. 1. — Figure illustrating how the PSF function J-(X,p' \p) may be used to accurately model a 
realistic mirror of finite extent. The top portion of the figure shows a "realistic" mirror with the 
detector in a de-focused position. In this position photons from a distant point source enter the 
mirror at the front and exit the mirror at the back, and are focused towards a point behind the 
detector causing the point source to appear out of focus. The bottom part of the figure shows that 
this effect can be modeled in terms of the concept of a "mirror node" , where rays enter and leave the 
same point with the appropriate redistribution J-(X,p' ,p) of "incoming" photons with direction p 
into "outgoing" ones with direction p' . In particular, note that the same PSF function appropriate 
for the de-focused position of the detector cannot be used to model the image when the detector 
is at the focal position; there a different PSF will be required. In other words, for F(X,p'.p) 
to accurately model the mirror, it also depends implicitly upon the location of the detector with 
respect to the focal surface. 
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assumed to contain the effects of (static) bad pixels, detector dead regions, detector boundaries, 
and so on, all of which cause it to vary spatially. 

Using the response function given by equation (5), the expected count rate from the source 
S(X,p) is seen to be 

S D (h,a,t) = T(a,t) J dxj dp R(h,a, X,p t ,t)S{X,p) (11) 
= T(a,t) J dX D(h, a, X) j dp' S(a - a(p' ,t)) j dp F(X,p' ,p t )M(X,p t )S(X,p), 

where, for notational simplicity, p t symbolizes TZ(t) ■ p. 

Telescope pointing motion will cause events to appear spatially mixed together when expressed 
as an image in detector coordinates. For this reason, it is preferable to work in the sky coordinate 
system where one can remove the effects of the motion by projecting the events to the sky in the 
appropriate manner. Thus, we define an aspect- corrected count rate via 

S A (h,p,t) = J(a(p t ,t),p t ,t)S D (h,a(p t ,t),t), (12) 

where J(a, p, t) is the instantaneous Jacobian of the transformation from the detector coordinate 
a to the sky coordinate p via the inversion of a(p t ,t). Physically, the Jacobian represents the 
stretching or magnification of an element of area on the detector as it appears in the sky. 

By making use of the well known change of variable formula for delta functions expressed in 
the form 

S(a(p t , t) - a(p>, t)) = J( [f'~ Pt } , (13) 
one can show that the aspect-corrected count rate is 

S A (h,p',t)=T(a(p' t ,t),t) j dX D(a(p' t , t), h, X) j dp f(X,p' t ,p t )M(X,p t )S(X,p). (14) 

The above equation assumes that complete knowledge of the aspect history is available in order to 
perform the aspect correction. In general, there will be uncertainties in the aspect solution which 
in turn leads to spatial uncertainties in the aspect corrected count rate. Mathematically, this will 
manifest itself as a broadening or smearing of the delta function in equation (13) by an amount 
that depends upon the aspect uncertainties. The most straightforward way to handle this effect is 
to absorb the uncertainties into the PSF itself. For this reason, in the following !Fa(X,p',p) will 
denote the PSF that includes the effect of the aspect uncertainties. 



3. Derivation of the ARF 



The total number of expected counts with pulse-height h over an observation interval r in 
some region can be computed by integrating t over the observation interval and pi over the sky 
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region f2, i.e., 

C n (h) = [ dp' [ T dtS A (h,p',t) (15) 
Jn Jo 

dX [ dp' [\tT(a(p' t ,t),t)D(a(p' t ,t),h,X)T A (X,p' t ,q t )M(X,q t )S(X,p). 
Jn Jo 

By using equation (8) and assuming for the moment that the RMF does not vary with position, 
the previous equation can be rewritten as 

Cn(h) = T cS J dXD R (h,X) J dp A n (X,p)S(X,p), (16) 

where 

A n (X,p) = — f dp' [ T dtT(t)Q(X,a(p' t ,t))F A (X,p' t ,p t )M(X,p t ) (17) 
T cB Jn Jo 

and the effective exposure time is given by 

r cff = f dtT(t). (18) 
Jo 

For simplicity it has been assumed that the good-time interval function T(t) does not depend upon 
detector position. 

Equation 17 defines the ARF. It has the units of area x counts per photon and depends upon 
the region fi, wavelength A, and sky position p. However, for the purposes of point-source analysis, 
knowledge of the ARF is required only at the position of the point source, where it may be regarded 
as depending only upon wavelength with the understanding that it is valid only for that source 
position and region. But for arbitrary sources, it should be regarded as an explicit function of p. 

Armed with the ARF, equation (16) is an integral equation that may be "solved" to yield the 
source distribution S(X,p) from the observed aspect-corrected counts Cn(h). Actually, because of 
the integration over the sky region $7, much of the spatial dependence in S(X,p) will be lost and 
in practice one will have to assume a known spatial dependence; the examples below illustrate this 
point more fully. One should also realize that the kernel of equation (16) is really a probability 
distribution and that the observed number of counts will most likely differ from the expected 
number of counts predicted by the equation. This implies that equation (16) does not really have 
a unique solution (for a finite observation time), and any method of "solving" should allow for 
fluctuations in the number of counts. One must also take into account any external background 
sources as well as any internal background produced by, say, noise in the detector. Techniques for 
treating this equation are beyond the scope of this paper and may be found elsewhere (Arnaud 
1996; Kahn and Blissett 1980). 

In the next few sections, equation (16) is considered in the context of various source distribu- 
tions. The problem of the practical computation of the ARF is taken up in the appendix. 
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3.1. A Single Point Source 

A point source located at position q in the sky with a spectrum s(A) may be represented using 
the source distribution 

S(X,p) = s(X)S(p-q). (19) 
Substituting this equation into equation (16) yields 

C Q (h) = T cS J dXD R (h,X)A n (X,q)s(X). (20) 

This equation is essentually the integral equation that the popular spectral analysis program xspec 
(Arnaud 1996) is designed to solve. As noted above, the ARF is required to be computed only 
at the position q of the point source; due to the time-dependence of the telescope motion, the 
integration over t in equation (17) does however sample various detector regions and variations in 
the off-axis effective area. 



3.2. Multiple Point Sources 

Multiple point sources may be represented by a source distribution of the form 



S(A,p) = J>(A)<5(p-<y, 



(21) 



where the position of the ith source is given by and its spectrum is Sj(A). Insertion of this 
distribution into equation (16) produces 



Ca{h) = r cB J dX D R (h, A) Yl M\ Qi)siW 



(22) 



If all of the sources have an identical source spectrum such that Sj(A) = s(A), then the resulting 
integral equation reduces to the case of a single point source, i.e., 



Cn(h) = r eS J dXD R (h,X) Y^AnMt 



s(X). 



(23) 



However, the more interesting case of spectrally distinct sources is more complicated to solve. In 
fact, its solution would require N integral equations since there are N unknown spectral distribu- 
tions Si(X). The most straightforward way to obtain the required number of independent equations 
would be to use N different regions Qi, not necessarily disjoint, and solve the resulting linear system 
of equations 

(24) 



C^{h)=r eS J dXD R (h,X) ^^(A,^-K(A) 



These equations are "coupled" to the extent that the Aq^X, qj) for i ^ j are non-zero, i.e., whether 
or not source j has a PSF contribution to region Qi (See figure 2). Such a system of equations may 
be handled using sherpa (Doe, et al. 1998), the Chandra Data System spectral analysis program. 
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Fig. 2. — Figure illustrating the use of multiple regions in equation (24) for a system of four point 
sources. In this figure, the thin radial lines emanating from each of the point sources represent the 
PSF (100% enclosed power) for the corresponding source, and the elliptical borders outline regions 
Qi over which counts are summed. The PSF from source s\ does not contribute to regions 

and Q4; hence, it effectively decouples from the other sources. Similarly, S4 may be treated by 
itself since none of the other sources contribute any counts to Q4. However, in dealing with S2, the 
contributions from both S3 and S4 must be taken into account. 



3.3. Extended Source with Uncorrelated Spatial and Spectral Distributions 



One of the simplest examples of an extended source is one in which the spatial and spectral 
distributions are uncorrelated. That is, the source distribution S(X,p) factors according to 

S(X,p) = s(X)p(p), (25) 

where p(p) defines the spatial distribution, assumed to be properly normalized such that 

1 = J dpp(p). (26) 

Combining this distribution with equation (16) yields 



Cn(h) = t cS J dA D R (h, A) J dp A Q (X,p)p(p) 



s(X) (27) 



Suppose that the form of p(p) is known. Then, the ARF could be combined with the known spatial 
distribution by defining 

A^(X) = J dpA n (X,p)p(p), (28) 
which leads to the xspec style equation 

Cn(h) = t cS j dA D R {h,X)A^(X)s(X) (29) 

for the unknown spectral function s(X). Of course, this methodology cannot be used if the spatial 
distribution p(p) is not known. The more general problem is addressed below in section 4 of this 
paper. 



3.4. An ARF in the presence of a spatially varying RMF 

It is important to note that the ARF, given in equation (17), is useful only when one can 
disregard spatial variations in the RMF. Unfortunately, this may not always be possible. For 
example, Chandra ACIS CCDs have a spatially varying response that must be properly taken into 
account. Provided that one wants to stay within the confines of the existing ARF+RMF paradigm, 
the only way to properly handle such cases is to filter the observed events over a region in detector 
coordinates where spatial variations in the RMF may be neglected. Mathematically, this procedure 
may be stated as follows. Let T denote the region on the detector where the RMF does not vary. 
Then define a filter F T (a) on this region by 



(30) 
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Multiplication of equation (11) by this filter, followed by aspect correction and summing over the 
sky region Q, yields 

Cl(h) = r cB j dA D T R (h, X) J dpA r n (X,p)S(X,p), (31) 

where 

A r n (X,p) = — f dp' [\tT(t)F r (a(p' t ,t))Q(X,a(p' t ,t))^ A (X,p' t ,pAM(X,p t ). (32) 
T cff Jo. Jo 

In these equations, C^(h) is the expected number of counts in the sky region f2 that also falls 
within the detector region T, A^(X,p) is the ARF appropriate for this region of the detector, and 
D^(h, A) is the region-dependent RMF. 



4. Definition of the Exposure Map 

The ARF presented in the previous section is primarily of use for spectral analysis over small 
spatial regions, e.g., the analysis of point sources. Much of the spatial information useful for the 
treatment of extended sources was lost in the construction of the ARF by integrating the response 
over a region of the sky. A related product, the exposure map, does not depend upon the 
integration over a sky region permitting it to be used for certain types of extended source analysis. 
The goal of this section is to define an exposure map and show how it may be used with extended 
sources. 

Start by integrating equation (14) over an observation time r to produce 

C(h,p') = f T dtS A (h,p',t) (33) 
Jo 

= £dtj dXT(a(p' t ,t),t)D(a(p',t),h,X) J dp F A (X,p' t ,p t )M(X,p t )S(X,p). 

Here, C(h,p') represents the expected total number of aspect-corrected counts with pulse-height h 
attributed to the sky position p' . 

In general the PSF is small for rays not too far off the optical axis, although it can become 
quite large for far off-axis rays. Suppose that the pointing motion amplitude is small enough that 
the PSF may be regarded as a scalar under the motion, i.e., 

F A (X,p',p)=F A (X,p' t ,p t ), (34) 

and then consider the integration over p in equation (33). For off-axis positions where the size of the 
PSF is small, for a fixed p', only a narrow range of p contributes to the integral. For the moment, 
assume that the mirror effective area does not vary much over this range. Then the approximation 

J dpF A (X,p' t ,p t )M(X,p t )S(X,p)^M(X,p>) J dpT A (X,p',p)S(X,p) (35) 
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can be used in equation (33) to yield 

C(h,p') = JdX dtT(a(p' t ,t),t)D(a(p>,t),h,\)M(\,p>)SA\p'), (36) 
where a PSF-smeared source Sy(A,p') has been defined by 

Sr(\,&) = j dpT A (X,p',p)S(X,p). (37) 
With the introduction of the exposure map E(h,X,p) defined by 

E(h,X,p) = — f dtT(a(p t ,t),t)D(a(p t ,t),h,X)M(X,p t ), (38) 

T eff JO 

equation (36) may be recast as 

C(h,p)=T cB J dXE(h,X,p)Sr(X,p). (39) 

It is important to understand that this is an integral equation describing the PSF-smeared source 
and not the true source. After "solving" this equation, one still has the task of removing the 
effects of the PSF to determine the true source. Nevertheless, equation (39) does have one very 
important feature not shared by the equations involving the ARF; namely, the spatial distribution 
of the expected aspect corrected counts, C(h,p), is the same as the PSF-smeared source's spatial 
distribution. 

A common use of the exposure map is to remove instrumental artifacts in images to obtain a 
better looking image. This is also known as "flux-correcting" the image. The method essentially 
assumes that the pulse-height resolution of the detector permits the separation of counts originating 
from photons from different energy bands, supplemented by the assumption that the source flux may 
be regarded as constant within a band (Snowden et al. 1994). To express this idea in quantitative 
terms, consider a range of Ah of pulse-heights centered on some pulse-height h and assume that the 
RMF A) is such that only those photons from the wavelength band A A about A can produce 

pulse-heights in the specified range. Now sum equation (39) over this range and consider only the 
photons from the wavelength band AA to yield 



C(h,p) = r off f dX Y, E(h,\,P)SAKP), ( 4 o) 
heAh JxeAX h£Ah 

which may be written in the more compact form 

C(Ah,p) = r eS f dXE(Ah,X,p)S r (X,p), (41) 

JXpAX 



'xeAX 

with the understanding that the pulse- height range Ah is summed over. 
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If the bandwidth AA is such that E(Ah, A, p) does not vary much over the band, then it may 
be removed from the integrand to obtain 



This equation says that the integrated profile of the PSF-smeared source flux over the wavelength 
band may be obtained by dividing the exposure map into the counts image constructed from the 
appropriate pulse-height range. 

The resolution of the source spectrum obtained by this technique is generally poorer than that 
of the detector's energy resolution because the wavelength band AA must be large enough to cover 
all the wavelengths that could contribute to the range of pulse-heights Ah. At the same time it 
must be small enough to ensure that E(Ah, X,p) may be treated as a constant in the wavelength 
band. It is possible that there may be bands in which these constraints are mutually exclusive. 
For this reason, the use of the exposure map is limited to situations where spectral resolution is 
of secondary importance. For example, spatial resolution is much more important than spectral 
resolution when doing source detection. For such a situation, one would run the source detection 
algorithm on an image obtained by dividing the total counts image over an exposure map integrated 
over the bandpass of the telescope. Another application of the exposure map would be to use it to 
get a crude estimate of the true spectrum, and use that as the first approximation in some more 
refined technique. 

Before leaving this section, it is important to point out that equation (39) is valid only as long 
as the size of the PSF is small enough that any variation in the effective area over the PSF can 
be neglected. If this is not the case, then it is impossible to give a definition of an exposure map 
that has the simple relationship between the observed counts image and the PSF-smeared source as 
described by this equation. By implication it follows that such an exposure map cannot be used for 
flux correction via the simple division of equation (42). However, if the detector response D(a, h, A) 
is uniform, then it is possible to commute the response with F(\,p' ,p) in equation (33) to produce 



This means that one must first deconvolve the effects of the PSF before correcting with the exposure 
map. The feasability of this will depend upon the energy dependence of the PSF where one may 
have to perform the deconvolution in specific energy bands. This prescription resembles the one 
advocated by White and Buote (2000) for the analysis of ASCA data. Alternatively, Ikebe (1995) 
has argued that one start essentually from equation (33) and employ a "forward-folding" method 
to estimate the source distribution S(X,p). 




(42) 




(43) 
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5. Derivation of the Grating ARF 

When diffraction gratings are added, the subsystem factorization, equation (5), must be mod- 
ified to 



R(h,a,X,p,t) = D(h,a,X) J dp' 5{a - a G (p' ,t)) ^5m(A)J ?r m (A,p / ,p) 



M(X,p), (44) 



where m represents the diffraction order and 8(a — cr G (p',t)) symbolizes the coordinate transfor- 
mation of a diffracted ray at the grating node with direction p' to the detector coordinate a. The 
actual diffraction into the rath order is represented by the term g m (X)J r m (X,p' ,p), which gives the 
probability for a ray with direction p and wavelength A to diffract into the mth order with direc- 
tion p'. The function g rn {\) is the mth order grating efficiency, and by definition the redistribution 
function satisfies the normalization condition 

1 = J dp' F m (X,p',p). (45) 

Despite the similarity in form of J-" m (X, p',p) to the imaging PSF F(X,p' ,p), it is important to 
appreciate one very important difference between these two functions. The imaging PSF is sharply 
peaked about the set of directions p' near p. However, J r m (X, p' ', p) is sharply peaked about a set 
of directions p' that vary linearly with wavelength A according to the diffraction equation 

(p'-p)xn = ^i, (46) 

where d is the grating period, the vector h is normal to the plane of the grating, and I is in the 
direction of the grating bars (see figure 3). 

It follows trivially from equation (44) that the expected count-rate into mth order is given by 



S { ™\h,a,t) =T(a,t) J dX 



D(h,a,X)g m (X) (47) 
x J dp' 5(a-a G (p',t)) J dpr m (X,p',p t )M(X,p t )S(X,p) . 

As in the imaging case, an aspect-corrected count-rate may be defined by 

S { ™\h,p,t) = J{a{p t ,t),p t ,t)S% n \h,a{p t ,t),t), (48) 

with the result 

S^\h,p',t)=T(a G (p' t ,t),t) J dXD(a G (p' t ,t),h,X)g m (X) j dp F m (X,p' t ,p t )M(X,p t )S(X,p). 

(49) 

This equation may be simplified for the special case of a point source located at q with a 
spectrum s(X), i.e., 

S(X,p') = s(X)5(p' -q). (50) 
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Fig. 3. — Figure illustrating how the grating RMF, Gq 1 \\, q), represents a redistribution from 
wavelength A to regions Q{. In this figure, the triad of unit vectors (n, I, d) specify an orthonormal 
coordinate basis centered upon the diffraction grating, with h normal to the surface of the grating 
and I in the direction of the grating bars. For a perfect grating, a photon with direction q and 
wavelength A will diffract into the direction p in accordance with equation (54) and a specified 
diffraction order m. For a realistic grating, the diffraction process must be described by a probability 
distribution T m (\,p, q) that is sharply peaked around the set of values that satisfy equation (54). 
The mth order grating RMF, Gq*\\, q) defined by equation (56), represents the redistribution 
probability for an incoming photon with direction q and wavelength A to diffract into a region 
In other words, for a specified zeroth order direction q, the grating RMF may be regarded as a 
mapping from wavelength A to region fij in much the same way as the detector RMF represents a 
mapping from wavelength A to pulse-height h. 
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In addition, assume that the telescope pointing motion amplitudes are small enough that the 
grating redistribution function behaves like a scalar under the motion as in equation (34). Then, 
the integral over pi may be readily performed to yield 

S^\h,p,t) = T(a G (p t ,t),t) J dXD(a G (p t ,t),h,X) 9rn (X)F m (X,p,q)M(X,q t )s(X). (51) 

Integration of this equation over a set of regions J2j and time r yields for the total expected number 
of counts with pulse-height h in the regions 

C&\h)= fdX f dpL(A) [ T dtT(a G (p t ,t),t)D(a G (p t ,t),h,X)M(X,q t )]^ m (X,p,q)s(X). 
J JQi L Jo 

(52) 

For fixed A and q, the sharp peaked nature of the grating redistribution function F m {X, p, q) implies 
that only a very narrow set of directions p will contribute to the term in square brackets. Moreover, 
any telescope pointing motion is expected to smooth out any non-uniformities in the detector QE 
appearing in this term such that one can evaluate it using the value of p determined by the grating 
equation. In other words, the term in square brackets can be replaced by a function A m (h, X) 
defined by 

A m (h, A) = g m (X)— [ dt T(a G (p t , t),t)D(a G (p t ,t), h, A)M(A, q t ), (53) 

T cS Jo 

where p satisfies 

(p-q)xh = ^i. (54) 
Hence, the number of counts in the mth order with pulse-height h is expected to be 

C^(h) = r cH J dX G^\x,q)A m (h,X)s(X), (55) 

where 

G^\x,q)= / dpF m (X,p,q), (56) 
jQi 

and T eff is given by equation (18). 

For reasons that will soon become clear, A m (h, A) is called the grating ARF, and Gq*\x, q) is 
called the grating RMF. To see this, consider the meaning of equation (56). For fixed A, equation 
(56) represents a redistribution from wavelength A to the region f2j, which may be regarded as 
the ith bin in p-space (see figure 3). This is the analog of the imaging RMF which describes a 
redistribution from A to a bin in pulse-height space. With this interpretation, equation (55) is 
formally identical the to equation equation (20), provided that one identifies A m (h, A) with the 
ARF. Hence, any techniques that are applicable to equation (20) may be readily applied to the 
solution of equation (55). 

Although A m (h, A) depends upon the pulse-height, in practice events will be filtered upon the 
pulse-height in order to perform order separation, provided that the intrinsic energy resolution of 
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the detector is adequate. For detectors with poor energy resolution, some other means of identifying 
mth order events will have to be used. In any case, A m (h, A) will most likely be summed over the 
range of pulse-heights appropriate to mth order events. In fact, the pulse-height range will generally 
vary with the wavelength such that the quantity 

hi (A) 

An(A)= A ™(h,\) (57) 

h=h (X) 

will actually be what is used in practice. For this reason, it is preferable to define the summed 
quantity A m (X) as the grating ARF. 

6. Conclusion 

In this paper, explicit expressions for the imaging ARF, grating ARF, and the exposure map 
were given in terms of the underlying instrumental responses that are consistent with the current 
use of the objects. These quantities were obtained from first principles by relating the expected 
detector count-rate to an incident photon source flux via the overall telescope response function 
suitably factored into individual instrumental responses. 

One of the complications in the derivation of these quantities concerned the proper treatment of 
time- varying effects due to telescope pointing motion, e.g., dither. At the same time, the assumed 
presence of motion about some nominal pointing allowed some important factorizations to take 
place that otherwise would have been suspect in regions containing detector boundaries or bad 
pixels. For this reason, purposely dithering an observation is recommended, provided, of course, 
that one can reconstruct the aspect history with sufficient accuracy. 

An added benefit of the first principles approach taken here is that it allows one to consider 
problems that cannot readily be handled by conventional means through the use of an ARF and 
an RMF. For example, as shown in section 3.3, it ARF may be applied to the analysis of an 
extended source provided one knows a priori that the source flux distribution factors into a known 
spatial component and an unknown spectral component. It is easy to find sources where such a 
factorization is not permissible; the supernova remanent, Cassiopeia A, is one. Another problem 
that does not appear to be treatable through standard the techniques is the analysis of an extended 
source in the presence of a diffraction grating. The grating ARF defined in section 5 was derived 
assuming a point source distribution. The basic problem with the analysis of an extended source 
is that, unlike a point source, there is no unique zeroth order position that one could use in the 
grating equation. By judicious filtering in pulse-height space, one may find regions where there is 
enough of a point-like behavior to permit the grating ARF to be used. However, how to handle a 
generic extended source in the presence of a diffraction grating is still an open question. It is hoped 
that the mathematical formulation of the extended source problem as given in sections 3 and 5 will 
lead to better insights into these problems and ultimately to their solution. 



This work also highlights some important practical considerations that should be taken into 
account in the design of astronomical data analysis software systems. For instance, to allow for 
the possibility of spatial variation in the underlying detector redistribution function, the software 
component responsible for the filtering of events should allow the user to easily filter simultaneously 
on both sky coordinates and detector coordinates. In addition, both filters would need to be passed 
to the program that generates the ARF. Finally, spectral fitting programs should be enhanced to 
facilitate the analysis of blended sources by handling the coupled integral equations in equation 
(24). 

I am especially grateful to David Huenemoerder for clarifying a number of issues during the 
course of this work. In addition, I also thank David Davis, John Houck, Norbert Schulz, and 
Michael Wise for useful discussions. Finally, I am indebted to Dan Dewey for his critical reading 
of the text and valuable suggestions regarding the presentation of some of the material. This work 
was supported under Chandra X-Ray Center contract SV1-61010 from the Smithsonian Institution. 



A. Numerical Considerations 

In this appendix, some "approximations" used for the practical computation of the ARF and 
grating ARF are discussed. In fact, these approximations are actually employed by the Chandra 
exposure map code suite for the generations of exposure maps and ARFs. Since the code is freely 
available 3 , the actual implementation details will not be discussed here. The reader should also 
note that some of these approximations may only be valid for the Chandra telescope, which dithers, 
and for other missions one may have to resort to the full definitions given in the main body of the 
text. 



A.l. Performing the Time Integrations via an Aspect Histogram 

The integrals over the observation time appearing in the equations for the ARF and the grating 
ARF can be quite computationally expensive, especially for long observation times. The general 
form of these integrals is given by 

1=1 dtT(t)F(a(t)) (Al) 
Jo 

where a(t) is an N dimensional time-dependent vector that characterizes the dither of the telescope 
and the relative motion of its subsystems. For example, A is 3 if there is no internal movement, 
and the dither is characterized by the roll, pitch, and yaw of the telescope. 



3 See http://chandra.harvard.edu for more information. 
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By multiplying the preceeding equation by the identity 

1 = daS(a- a(t)) 



(A2) 




(A3) 



where 



U{a) = ! dt T(t)5(a - a(t)) 



(A4) 



The quantity TC(a)da has a very simple interpretation. It represents the total amount of time, 
weighted by T(t), that the point spent in the volume element Act at ct. Now, if the telescope 
dithers around some nominal pointing, and if the time-dependent internal motions due to, e.g., 
thermal expansion are small, then the point ot{t) will be confined to some small volume in the N 
dimensional space. This means that Tt(a) will be non-zero only in that small volume and zero 
everywhere else. So, to compute the time integration over long observation times for the case 
of small dither amplitudes, it is often more efficient to compute the value of H(a) and use it to 
evaluate equation (A3). In practice, the portion of the N dimensional space where H(a) is non-zero 
is sub-divided into small volume elements Act. Then the discretized quantity Ha-Aa = H(at)Aat 
is computed and used in a discretized version of equation (A3), i.e., 



For reasons that should be apparent, Ha-Aa is called the aspect histogram. 

The major advantage of this approach for the case of small dither amplitudes is that there are 
likely to be many fewer terms to sum in equation (A5) than if a straightforward discretization were 
used to perform the time integration in equation (Al). Moreover, there are efficient algorithms 
based upon 2^-ary trees for computing the aspect histogram. For example, the code for both the 
Chandra and ROSAT missions use an octtree for N = 3. (For Chandra, the value of N used is 3 
rather than 6 through the use of "effective" offsets.) 



The ARF is a complicated function requiring complete knowledge of the detector QE, mirror 
effective area, aspect solution, and the point spread function. To compute it directly from equation 
(17) or from equation (32) in the case of a spatially varying RMF, one would need to carry out an 
integration over time as well as a 2-d integration over the sky region, and do this for every point 
in the sky. Clearly, this is not practical and in view of the fact that there will be uncertainties in 
the instrumental responses at this level of detail, such a calculation is unwarranted. Instead, one 
can make several simplifying approximations that permit the ARF to be computed in an economic 
manner. 




(A5) 



ct 



A. 2. Computation of the Imaging ARF 
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As written, equation (17) is valid for any motion of the spacecraft, including slew. However, 
here it shall be assumed that one is dithering about some mean pointing and that the scale of the 
dither is small enough that any variations in the PSF and the mirror effective area on this scale 
can be neglected. Therefore, equation (17) will be approximated by 

An(A,p) « — / dp'F A (\, (p'), (p»M(A, <p» [ T dt T(t)Q(X,a(p' t ,t)), (A6) 

Toff JQ JO 

where (p) represents the time-average of p t , i.e., 

(p) = — [ T dtT(t)K(t)-p. (A7) 

r off JO 

Similarly, define 

<Q(A,p')) = — dtT(t)Q(\,*(p' t ,t)) (A8) 

r off JO 

to be the time-averaged value of the QE. Then one can write 

An(A,p) « M(A, (p)) / dp'^ A (A, <p'), (p))(Q(A,p')). (A9) 
jn 

The time-averaging over the dither motion has the effect of smoothing out any large variations in 
the QE over the region. In fact, this is the primary purpose of the dither. Now since (Q(X,p')) can 
be assumed to vary slowly over the region, and since Fa(X,p' ,p) is expected to rapidly go to zero 
as p' moves away from p, (Q(X,p')) can be replaced by its average over the region and removed 
from the integrand. This leads to the result 

An(X,p) « /n(A,p)M(A, (p»(Q(A)) n , (A10) 

where 

(Q(A)}n = ^ / dp / (Q(A,p / )) (All) 
is the average of (Q(X,p)) over the region and the PSF fraction in the region is given by 

/n(A,p) = / dp'^ 4 (A,(p , ),(p)) (A12) 
Jn 

dp F A (\,p',p). 



/• 

Jn 



A. 3. Computation of the Grating ARF 

The grating ARF is defined by equation equation (57), rewritten here as 

fci(A) 

A m (X) 



J2 D R (h,X) xg m (X)— [ dtT(t)Q(a G (p t ,t),t)M(X,q t ), (A13) 

_, J r off JO 



h=h (\) 
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where equation (8) has been used with the assumption that the RMF does not vary spatially. For 
a spatially varying RMF an additional spatial filter would need to be applied as was done for the 
imaging ARF to derive equation (32). In the above equation, p depends upon the source position 
and the wavelength according to 

(p-q)xh = —I. (A14) 

If the amplitude of the dither is small on the scale of the variations in the mirror effective area, 
then M(A, q) may be replaced by M(A, (q)) and removed from the integrand. Hence, the grating 
ARF may be approximated by 

hi (A) 

A m (X)^g m (X)(Q(X,p))M(X,(q}) £ D R (h, A), (A15) 

h=h (X) 

where 

(Q(A, p)> = — f T dt T(t)Q(a G (p t ,t),t). (A16) 

T cS Jo 
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